This document contains a description of optimisation tools for electric system operation simulation. If you are interested in planing you might prefer to switch directly to file optim-Planing.ipynb. In files BasicFunctionalities/input-XXX.ipynb you can learn to understand input data (consumption, availability).
This document will gives a chance to understand
It proposes to enter the subject by increasing progressively the number of variables and constraints in the optimisation problem, hence moving toward more realism through the document, introducing:
It relies on different test cases that allow to
If, after reading this file, you want to build your own pyomo model you can go to optim-Operation-Advanced.ipynb.
#region importation of modules
import os
if os.path.basename(os.getcwd())=="BasicFunctionalities":
os.chdir('..') ## to work at project root like in any IDE
InputFolder='Data/input/'
import numpy as np
import pandas as pd
import csv
import datetime
import copy
import plotly.graph_objects as go
import matplotlib.pyplot as plt
from sklearn import linear_model
import sys
from functions.f_operationModels import *
from functions.f_optimization import *
from functions.f_graphicalTools import *
#endregion
#region Solver and data location definition
InputFolder='Data/input/'
if sys.platform != 'win32':
myhost = os.uname()[1]
else : myhost = ""
if (myhost=="jupyter-sop"):
## for https://jupyter-sop.mines-paristech.fr/ users, you need to
# (1) run the following to loanch the license server
if (os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmgrd -c /opt/mosek/9.2/tools/platform/linux64x86/bin/mosek.lic -l lmgrd.log")==0):
os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmutil lmstat -c 27007@127.0.0.1 -a")
# (2) definition of license
os.environ["MOSEKLM_LICENSE_FILE"] = '@jupyter-sop'
BaseSolverPath='/Users/robin.girard/Documents/Code/Packages/solvers/ampl_macosx64' ### change this to the folder with knitro ampl ...
## in order to obtain more solver see see https://ampl.com/products/solvers/open-source/
## for eduction this site provides also several professional solvers, that are more efficient than e.g. cbc
sys.path.append(BaseSolverPath)
solvers= ['gurobi','knitro','cbc'] # try 'glpk', 'cplex'
solverpath= {}
for solver in solvers : solverpath[solver]=BaseSolverPath+'/'+solver
solver= 'mosek' ## no need for solverpath with mosek.
#endregion
In this section we propose a simple version of the optimisation problem with a mathematical description. This allows us to discuss the use of pyomo (a python package to write optimisation problems).
Mathematical formulation as a linear programming problem
\begin{align} &\text{Cost function }& &\min_{x} \sum_t \sum_i \pi_i x_{it}\;\;\; & & \pi_i \text{ marginal cost}\\ &\text{Power limit } & &\text{ s.t.} \;\; 0 \leq x_{it}\leq a_{it} \bar{x_i} & &\bar{x_i} \text{ installed power, } a_{it} \text{ availability}\\ &\text{Meet demand } & & \sum_i x_{it} \geq C_t && C_t \text{ Consumption}\\ &\text{Stock limit } & &\sum_t x_{it}\leq E_i && E_i=\bar{x_i}*N_i \text{ Energy capacity limit}\\ \end{align}This is linear programing and could be transformed into something like
$ \begin{align} & \min_y & c^Ty \\ & Ay\leq b \end{align} $
with a well chosen parameter matrix $A$, parameter vector $b$ and $c$ and variable vector $y$. While increasing complexity of the problem, this can become very painful.
With the Python package Pyomo it is almost sufficient to write the mathematical equations. Pyomo is then charged of building the matrix form, and you just have to think about problem formulation. Pyomo is similar to other tools such as GAMS, AMPL, ...
Principles of a pyomo model To build a model, Pyomo needs three different kinds of data : sets, parameters and variables.
Sets are dimensions, here the time and the name of technology plus a mix of these two : TIMESTAMP, TECHNOLOGIES and TIMESTAMP_TECHNOLOGIES (product set).
Parameters are tables indexed by set whose values are specified by the user. Here is the list of the parameters we use in the simplest cases : energycost, EnergyNbhourCap, capacity, availability factor, area consumption.
Variable are tables indexed by set whose values are found by the solver : the energy produced by each mean of production.
Now if you wish to learn more about pyomo and to see how the optimisation problem is writen in pyomo language, take look at function GetElectricSystemModel_GestionSingleNode (by doing a control-click on it in the code that follos a bit later, otherwise it is in functions/f_operationModels.py). We will use this function in the following exemple.
First simulation -- loading parameters
We start with the simple test case bellow, with only two production means (nuke and thermal) but you can run this with more than that, and change the installed power. If you want to change assumptions you will need to change the csv files that are loaded below.
#region I - Simple single area : loading parameters
Zones="FR" ; year=2013
#### reading areaConsumption availabilityFactor and TechParameters CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Gestion-Simple_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0).set_index(["TECHNOLOGIES"])
#### Selection of subset
Selected_TECHNOLOGIES=['OldNuke','CCG'] #you can add technologies here
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
TechParameters.loc["CCG",'capacity']=100000 ## margin to make everything work
#TechParameters.loc["WindOnShore",'capacity']=117000
#TechParameters.loc["Solar",'capacity']=67000
#endregion
#endregion
First simulation -- running the optimisation
Now we run the optimisation script and load optimised variables.
At this stage, or later, you might be willing to to see the pyomo model. If it is the case, dig into function GetElectricSystemModel_GestionSingleNode (by doing a control-click on it, otherwise it is in functions/f_operationModels.py)
#region I - Simple single area : Solving and loading results
model = GetElectricSystemModel_GestionSingleNode(areaConsumption,availabilityFactor,TechParameters)
if solver in solverpath : opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
## result analysis
Variables=getVariables_panda_indexed(model)
#pour avoir la production en KWh de chaque moyen de prod chaque heure
production_df=Variables['energy'].pivot(index="TIMESTAMP",columns='TECHNOLOGIES', values='energy')
### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
print(production_df.sum(axis=0)/10**6) ### energies produites TWh
print(Variables['energyCosts']) #pour avoir le coût de chaque moyen de prod à l'année
#endregion
TECHNOLOGIES CCG 80.112314 OldNuke 411.902306 dtype: float64 TECHNOLOGIES energyCosts 0 CCG 5.511727e+09 1 OldNuke 2.883316e+09
We have analysed the result of the optimisation first with the dictionary variables that contains all optimized variables. The most important variable here is 'energy'. The other one 'energyCosts' can be computed from 'energy' and it is somehow redundant.
Lagrange multiplier might be more difficult to understand for those who still lack a good optimisation course (the one at second semester of first year of MINES ParisTech is perfect). In cases you want to dig this, you can have a look at Boyd's course, e.g. starting with lecture 8 or before. The most important message here is that lagrange multipliers associated to the demand constraint (here called 'energyCtr') are meant to mimic market prices. Lagrange multiplier associated to this constraint at time t is the marginal cost that one would pay to increase $C_t$ by a small amount. They can be used to simulate market prices.
#region I - Simple single area : visualisation and lagrange multipliers
### representation des résultats
TIMESTAMP_d=pd.date_range(start=str(year)+"-01-01 00:00:00",end=str(year)+"-12-31 23:00:00", freq="1H")
production_df.index=TIMESTAMP_d; areaConsumption.index=TIMESTAMP_d;
fig=MyStackedPlotly(y_df=production_df,Conso = areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#### lagrange multipliers
Constraints= getConstraintsDual_panda(model)
# Analyse energyCtr
energyCtrDual=Constraints['energyCtr']; energyCtrDual['energyCtr']=energyCtrDual['energyCtr']
energyCtrDual
round(energyCtrDual.energyCtr,2).unique()
# Analyse CapacityCtr
CapacityCtrDual=Constraints['CapacityCtr'].pivot(index="TIMESTAMP",columns='TECHNOLOGIES', values='CapacityCtr');
round(CapacityCtrDual,2)
round(CapacityCtrDual.OldNuke,2).unique() ## if you increase by Delta the installed capacity of nuke you decrease by xxx the cost when nuke is not sufficient
round(CapacityCtrDual.CCG,2).unique() ## increasing the capacity of CCG as no effect on prices
#endregion
array([-0.])
We have the Lagrange multipliers associated to each (active) constraint (zero means unactive constraint) in dictionary 'Constraints'. We have three different constraints : the energy constraint, the capacity constraint and the storage constraint (unactive here) Try understanding the meaning of lagrange multipliers, here and in the next Sections.
In the this section, we will increase the complexity of the problem given in Section 2 and add : dependency on area z (country), a congestion constraint, ramp constraints.
\begin{align} &\text{Cost function }& &\min_{x} \sum_z \sum_t \sum_i \pi_{iz} x_{itz}\;\;\; & & \pi_{iz} \text{ marginal cost}\\ &\text{Power limit } & &\text{ s.t.} \;\; 0 \leq x_{itz}\leq a_{itz} \bar{x_{iz}} & &\bar{x_{iz}} \text{ installed power, } a_{itz} \text{ availability}\\ &\text{Meet demand } & & \sum_i x_{itz} \geq C_{tz} && C_{tz} \text{ Consumption}\\ &\text{Stock limit } & &\sum_t x_{it}\leq E_i && E_i=\bar{x_i}*N_i \text{ Energy capacity limit}\\ &\text{ramp limit } & &rc^-_i *x_{it}\leq x_{it}-x_{i(t+1)}\leq rc^+_i *x_{it} && rc^+_i rc^-_i\text{ ramp limit}\\ \end{align}#region II - Ramp Ctrs Single area : loading parameters loading parameterscase with ramp constraints
Zones="FR"
year=2013
Selected_TECHNOLOGIES=['OldNuke', 'CCG',"curtailment"] #you'll add 'Solar' after
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Gestion-RAMP1_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0).set_index(["TECHNOLOGIES"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
TechParameters.loc["CCG",'capacity']=100000 ## margin to make everything work
TechParameters.loc["OldNuke",'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc["OldNuke",'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#TechParameters.loc["WindOnShore","energyCost"]=0.001 ## a bit strong to put in light the effect
#endregion
#region II - Ramp Ctrs Single area : solving and loading results
model = GetElectricSystemModel_GestionSingleNode(areaConsumption,availabilityFactor,TechParameters)
opt = SolverFactory(solver)
results=opt.solve(model)
Variables=getVariables_panda_indexed(model)
#pour avoir la production en KWh de chaque moyen de prod chaque heure
production_df=Variables['energy'].pivot(index="TIMESTAMP",columns='TECHNOLOGIES', values='energy')
### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
print(production_df.sum(axis=0)/10**6) ### energies produites TWh
print(Variables['energyCosts']) #pour avoir le coût de chaque moyen de prod à l'année
#endregion
TECHNOLOGIES CCG 86.319599 OldNuke 405.657620 curtailment 0.037401 dtype: float64 TECHNOLOGIES energyCosts 0 CCG 5.938788e+09 1 curtailment 1.122030e+08 2 OldNuke 2.839603e+09
#region II - Ramp Ctrs Single area : visualisation and lagrange multipliers
### representation des résultats
TIMESTAMP_d=pd.date_range(start=str(year)+"-01-01 00:00:00",end=str(year)+"-12-31 23:00:00", freq="1H")
production_df.index=TIMESTAMP_d; areaConsumption.index=TIMESTAMP_d;
fig=MyStackedPlotly(y_df=production_df,Conso = areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
#### lagrange multipliers
Constraints= getConstraintsDual_panda(model)
# Analyse energyCtr
energyCtrDual=Constraints['energyCtr']; energyCtrDual['energyCtr']=energyCtrDual['energyCtr']
energyCtrDual
round(energyCtrDual.energyCtr,2).unique()
# Analyse CapacityCtr
CapacityCtrDual=Constraints['CapacityCtr'].pivot(index="TIMESTAMP",columns='TECHNOLOGIES', values='CapacityCtr');
round(CapacityCtrDual,2)
round(CapacityCtrDual.OldNuke,2).unique() ## if you increase by Delta the installed capacity of nuke you decrease by xxx the cost when nuke is not sufficient
round(CapacityCtrDual.CCG,2).unique() ## increasing the capacity of CCG as no effect on prices
#endregion
array([-0.])
Again you can have to look at lagrange multipliers. Try adding renewable production.
Try modifying the ramp constraint to generate negative lagrange multipliers.
Math here are in 3.1
#region III - Ramp Ctrs multiple area : loading parameters
Zones="FR_DE_GB_ES"
year=2016
Selected_AREAS=["FR","DE"]
Selected_TECHNOLOGIES=['OldNuke', 'CCG','WindOnShore',"curtailment"] #you'll add 'Solar' after #'NewNuke', 'HydroRiver', 'HydroReservoir','WindOnShore', 'WindOffShore', 'Solar', 'Curtailement'}
#### reading CSV files
TechParameters = pd.read_csv(InputFolder+'Gestion_MultiNode_DE-FR_AREAS_TECHNOLOGIES.csv',sep=',',decimal='.',comment="#",skiprows=0).set_index(["AREAS","TECHNOLOGIES"])
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["AREAS","TIMESTAMP"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["AREAS","TIMESTAMP","TECHNOLOGIES"])
ExchangeParameters = pd.read_csv(InputFolder+'Hypothese_DE-FR_AREAS_AREAS.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["AREAS","AREAS.1"])
#ExchangeParameters.loc[ExchangeParameters.AREAS=="FR",'maxExchangeCapacity']=90000 ## margin to make everything work
#ExchangeParameters.loc[ExchangeParameters.AREAS=="DE",'maxExchangeCapacity']=90000 ## margin to make everything work
#### Selection of subset
TechParameters=TechParameters.loc[(Selected_AREAS,Selected_TECHNOLOGIES),:]
areaConsumption=areaConsumption.loc[(Selected_AREAS,slice(None)),:]
availabilityFactor=availabilityFactor.loc[(Selected_AREAS,slice(None),Selected_TECHNOLOGIES),:]
TechParameters.loc[(slice(None),'CCG'),'capacity']=100000 ## margin to make everything work
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion
#region III - Ramp Ctrs multiple area : solving and loading results
### small data cleaning
availabilityFactor.availabilityFactor[availabilityFactor.availabilityFactor>1]=1
model = GetElectricSystemModel_GestionMultiNode(areaConsumption,availabilityFactor,TechParameters,ExchangeParameters)
opt = SolverFactory(solver)
results=opt.solve(model)
Variables=getVariables_panda_indexed(model)
production_df=EnergyAndExchange2Prod(Variables)
### Check sum Prod = Consumption
Delta= production_df.sum(axis=1)-areaConsumption.areaConsumption
abs(Delta).sum()
## adding dates
production_df_=ChangeTIMESTAMP2Dates(production_df,year)
areaConsumption_=ChangeTIMESTAMP2Dates(areaConsumption,year)
fig=MyAreaStackedPlot(production_df_,Conso=areaConsumption_)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df.groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df.groupby(by=["AREAS"]).max()
Constraints= getConstraintsDual_panda(model)
Constraints.keys()
Constraints['energyCtr']
#endregion
| AREAS | TIMESTAMP | energyCtr | |
|---|---|---|---|
| 0 | FR | 1 | 68.8 |
| 1 | FR | 2 | 68.8 |
| 2 | FR | 3 | 68.8 |
| 3 | FR | 4 | 68.8 |
| 4 | FR | 5 | 68.8 |
| ... | ... | ... | ... |
| 17563 | DE | 8780 | 68.8 |
| 17564 | DE | 8781 | 68.8 |
| 17565 | DE | 8782 | 68.8 |
| 17566 | DE | 8783 | 68.8 |
| 17567 | DE | 8784 | 68.8 |
17568 rows × 3 columns
Just have a look at optim-Storage.ipynb
This section is a bit more difficult that the other one. You can skip the math if you're only interested in using the storage optimisation tool.
Let us define $\psi : z\in \mathbb{R}^T \rightarrow \psi(z)\in \mathbb{R}$ as the solution to the simple operation problem
\begin{align} &\text{Cost function }& &\min_{x} \sum_t \sum_i \pi_i x_{it}\;\;\; & & \pi_i \text{ marginal cost}\\ &\text{Power limit } & &\text{ s.t.} \;\; 0 \leq x_{it}\leq a_{it} \bar{x_i} & &\bar{x_i} \text{ installed power, } a_{it} \text{ availability}\\ &\text{Meet demand } & & \sum_i x_{it} \geq C_t+z_t && C_t \text{ Consumption}\\ &\text{Stock limit } & &\sum_t x_{it}\leq E_i && E_i=\bar{x_i}*N_i \text{ Energy capacity limit}\\ \end{align}Adding one storage to this operation problem is done here : \begin{align} &\text{Cost function }& &\min_{z} \psi(z) \\ &\text{storage power constraint}& &-p_{\max}\leq z_i\leq p_{\max} \\ &\text{storage energy constraint}& &-p_{\max}\leq \sum_{j=1}^i z_j\leq p_{\max} \\ \end{align}
Note that $\psi(z)=\sum_i \psi_i(z_i)$. When we solve the planing problem with a value of storage operation $z^0$ this not only provides $\psi(z^0)$ but also the gradient $\nabla \psi(z^0)$ (with Lagrange multipliers), and if we denote $\delta^0_i=(\nabla \psi(z^0))_i$ it entails a convenient uniform lower bound on $\psi(z)$ :
$ \begin{align} &\beta^0(z)& = &\;\; \langle \nabla \psi(z^0),z-z^0\rangle +\psi(z^0) \\ & & = & \;\;\sum_i \;\;\delta^0_i z_i+ \psi_i(z_i^0) -\delta^0_iz_i^0 \\ & & \hat{=} & \;\;\sum_i \beta_i^0(z_i) \end{align} $
As a consequence, an iterative procedure is implemented to solve this problem,
Initialisation, Take
step $k=1...$
This algorithm is implemented in function GetElectricSystemModel_GestionSingleNode_with1Storage in file functions/f_operationModels.py. The idea of this algorithm is quite general and a version exists also in the multi-node case and for planing problems (in file functions/f_planingModels.py)
#region IV Ramp+Storage single area : loading parameters
Zones="FR"
year=2013
Selected_TECHNOLOGIES=['OldNuke','WindOnShore', 'CCG',"curtailment"]
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Gestion-RAMP1_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0).set_index(["TECHNOLOGIES"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
TechParameters.loc["CCG",'capacity']=100000 ## margin to make everything work
TechParameters.loc["OldNuke",'RampConstraintMoins']=0.02 ## a bit strong to put in light the effect
TechParameters.loc["OldNuke",'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
p_max=5000
StorageParameters={"p_max" : p_max , "c_max": p_max*10,"efficiency_in": 0.9,"efficiency_out" : 0.9}
#endregion
#region IV Ramp+Storage single area : solving and loading results
res= GetElectricSystemModel_GestionSingleNode_with1Storage(areaConsumption,availabilityFactor,
TechParameters,StorageParameters)
Variables = getVariables_panda_indexed(res['model'])
Constraints = getConstraintsDual_panda(res['model'])
areaConsumption = res["areaConsumption"]
production_df=Variables['energy'].pivot(index="TIMESTAMP",columns='TECHNOLOGIES', values='energy')
Delta= production_df.sum(axis=1)-areaConsumption["NewConsumption"]
sum(abs(Delta))
production_df.loc[:,'Storage'] = -areaConsumption["Storage"] ### put storage in the production time series
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df[production_df>0].sum(axis=0)/10**6 ### energies produites TWh
production_df.max(axis=0)/1000 ### Pmax en GW
TIMESTAMP_d=pd.date_range(start=str(year)+"-01-01 00:00:00",end=str(year)+"-12-31 23:00:00", freq="1H")
production_df.index=TIMESTAMP_d; areaConsumption.index=TIMESTAMP_d;
fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
stats=res["stats"]
#endregion
0 1 2
If you want to contribute here, try improving the visualisation tools (with respect to storage, multizone, integration of interconnexions ...)
#region V Ramp+Storage Multi area : loading parameters
Zones="FR_DE_GB_ES"
year=2016
Selected_AREAS=["FR","DE"]
Selected_TECHNOLOGIES=['OldNuke', 'CCG','WindOnShore',"curtailment"] #you'll add 'Solar' after #'NewNuke', 'HydroRiver', 'HydroReservoir','WindOnShore', 'WindOffShore', 'Solar', 'Curtailement'}
#### reading CSV files
TechParameters = pd.read_csv(InputFolder+'Gestion_MultiNode_DE-FR_AREAS_TECHNOLOGIES.csv',sep=',',decimal='.',comment="#",skiprows=0).set_index(["AREAS","TECHNOLOGIES"])
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["AREAS","TIMESTAMP"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["AREAS","TIMESTAMP","TECHNOLOGIES"])
ExchangeParameters = pd.read_csv(InputFolder+'Hypothese_DE-FR_AREAS_AREAS.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["AREAS","AREAS.1"])
#ExchangeParameters.loc[ExchangeParameters.AREAS=="FR",'maxExchangeCapacity']=90000 ## margin to make everything work
#ExchangeParameters.loc[ExchangeParameters.AREAS=="DE",'maxExchangeCapacity']=90000 ## margin to make everything work
#### Selection of subset
TechParameters=TechParameters.loc[(Selected_AREAS,Selected_TECHNOLOGIES),:]
areaConsumption=areaConsumption.loc[(Selected_AREAS,slice(None)),:]
availabilityFactor=availabilityFactor.loc[(Selected_AREAS,slice(None),Selected_TECHNOLOGIES),:]
TechParameters.loc[(slice(None),'CCG'),'capacity']=100000 ## margin to make everything work
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
p_max=10000
StorageParameters=pd.DataFrame([])
for AREA in Selected_AREAS :
StorageParameters_ = {"AREA": AREA, "p_max": p_max, "c_max": p_max * 10, "efficiency_in": 0.9,
"efficiency_out": 0.9}
StorageParameters=StorageParameters.append(pd.DataFrame([StorageParameters_]))
StorageParameters=StorageParameters.set_index("AREA")
#endregion
#region V Ramp+Storage multi area : solving and loading results
res= GetElectricSystemModel_GestionMultiNode_with1Storage(areaConsumption,availabilityFactor,
TechParameters,ExchangeParameters,StorageParameters)
Variables = getVariables_panda(res['model'])
production_df=EnergyAndExchange2Prod(Variables)
areaConsumption = res["areaConsumption"]
### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.NewConsumption); ## comparaison à la conso incluant le stockage
abs(Delta).max()
production_df.loc[:,'Storage'] = -areaConsumption["Storage"] #### ajout du stockage comme production
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()
production_df_=ChangeTIMESTAMP2Dates(production_df,year)
areaConsumption_=ChangeTIMESTAMP2Dates(areaConsumption,year)
fig=MyAreaStackedPlot(production_df_,Conso=areaConsumption_)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
abs(areaConsumption["Storage"]).groupby(by="AREAS").sum() ## stockage
production_df.groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df[production_df>0].groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df.groupby(by="AREAS").max()/1000 ### Pmax en GW ### le stockage ne fait rien en Allemagne ??? bizarre
production_df.groupby(by="AREAS").min()/1000 ### Pmax en GW
#endregion
0 1 2
| CCG | OldNuke | WindOnShore | curtailment | DE | FR | Storage | |
|---|---|---|---|---|---|---|---|
| AREAS | |||||||
| DE | 0.000000e+00 | 4.310174 | -0.0 | 0.0 | 0.0 | -9.0 | -9.0 |
| FR | -1.335820e-15 | 33.523576 | -0.0 | 0.0 | -9.0 | 0.0 | -9.0 |
If you want to contribute here :
#region VI Complete "simple" France loading parameters
Zones="FR"
year=2013
Selected_TECHNOLOGIES=['OldNuke','Coal','CCG','TAC', 'WindOnShore','HydroReservoir','HydroRiver','Solar','curtailment']
#### reading CSV files
areaConsumption = pd.read_csv(InputFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP"])
availabilityFactor = pd.read_csv(InputFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
sep=',',decimal='.',skiprows=0).set_index(["TIMESTAMP","TECHNOLOGIES"])
TechParameters = pd.read_csv(InputFolder+'Gestion-Simple_TECHNOLOGIES.csv',sep=',',decimal='.',skiprows=0).set_index(["TECHNOLOGIES"])
#### Selection of subset
availabilityFactor=availabilityFactor.loc[(slice(None),Selected_TECHNOLOGIES),:]
TechParameters=TechParameters.loc[Selected_TECHNOLOGIES,:]
#TechParameters.loc[TechParameters.TECHNOLOGIES=="CCG",'capacity']=15000 ## margin to make everything work
p_max=5000
StorageParameters={"p_max" : p_max , "c_max": p_max*10,"efficiency_in": 0.9,"efficiency_out" : 0.9}
#endregion
#region VI Complete "simple" France : solving and loading results
res= GetElectricSystemModel_GestionSingleNode_with1Storage(areaConsumption,availabilityFactor,
TechParameters,StorageParameters)
Variables = getVariables_panda_indexed(res['model'])
Constraints = getConstraintsDual_panda(res['model'])
areaConsumption = res["areaConsumption"]
production_df=Variables['energy'].pivot(index="TIMESTAMP",columns='TECHNOLOGIES', values='energy')
Delta= production_df.sum(axis=1)-areaConsumption["NewConsumption"]
sum(abs(Delta))
production_df.loc[:,'Storage'] = -areaConsumption["Storage"] ### put storage in the production time series
production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df[production_df>0].sum(axis=0)/10**6 ### energies produites TWh
production_df.max(axis=0)/1000 ### Pmax en GW
TIMESTAMP_d=pd.date_range(start=str(year)+"-01-01 00:00:00",end=str(year)+"-12-31 23:00:00", freq="1H")
production_df.index=TIMESTAMP_d; areaConsumption.index=TIMESTAMP_d;
fig=MyStackedPlotly(y_df=production_df, Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename='file.html') ## offline
fig.show()
stats=res["stats"]
#endregion
0 1 2